
!pip install beautifulsoup4
!pip install geojson
!pip install geopandas
!pip install keplergl
!pip install plotly
!pip install polyline
!pip install pyarrow
!pip install shapely
!pip install seaborn
!pip install scipy
!pip install tqdm


BASE_DIR = './'


from bs4 import BeautifulSoup
from datetime import time
from geojson import LineString, Feature, FeatureCollection
import geopandas as gpd
import json
import keplergl
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import pickle
import plotly.express as px
import polyline
import requests
from shapely import geometry as shapely_geometry
from tqdm import tqdm

px.set_mapbox_access_token('pk.eyJ1Ijoicm1vc3RlcnQiLCJhIjoiY2tiZ3JjdWpyMG9pOTJ5bXptazQxYzY3dyJ9.UxiB2rAg3oBNH0yx0L9LcA')


DATA_DIR = BASE_DIR + 'TfL Cycling Journey Data/'
INTERIM_OUTPUT_DIR = BASE_DIR + 'Intermediate Output/'
REPORT_OUTPUT_DIR = BASE_DIR + 'Report Output/'

for DIRECTORY in [DATA_DIR, INTERIM_OUTPUT_DIR, REPORT_OUTPUT_DIR]:
    if not os.path.exists(DIRECTORY):
        os.mkdir(DIRECTORY)


with open(DATA_DIR + 'cycling.data.tfl.gov.uk.html', 'r') as f:
    cycling_data_page = f.read()


soup = BeautifulSoup(cycling_data_page, 'html.parser')

# Pull out every Journey Data Extract file listed on the page:
journey_data_urls = []
for link_elements in soup.find_all('a'):
    link = (link_elements.get('href'))
    if 'JourneyDataExtract' in link:
        journey_data_urls.append(link)


journey_data_urls[:5]

['https://cycling.data.tfl.gov.uk/usage-stats/01aJourneyDataExtract10Jan16-23Jan16.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/01bJourneyDataExtract24Jan16-06Feb16.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/02aJourneyDataExtract07Fe16-20Feb2016.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/02bJourneyDataExtract21Feb16-05Mar2016.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/03JourneyDataExtract06Mar2016-31Mar2016.csv']


[url for url in journey_data_urls if '2022' in url][:5]

['https://cycling.data.tfl.gov.uk/usage-stats/298JourneyDataExtract29Dec2021-04Jan2022.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/299JourneyDataExtract05Jan2022-11Jan2022.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/300JourneyDataExtract12Jan2022-18Jan2022.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/301JourneyDataExtract19Jan2022-25Jan2022.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/302JourneyDataExtract26Jan2022-01Feb2022.csv']


SELECTED_YEARS = [2022]  # Filter years to be one of the years in this list
SELECTED_MONTHS = ['jun', 'jul']  # Filter months to be in this list

selected_journey_data = [
    url for url in journey_data_urls 
    if any([str(selected_year) in url for selected_year in SELECTED_YEARS]) 
    and any([selected_month.lower() in url.lower() for selected_month in SELECTED_MONTHS])
]


selected_journey_data

['https://cycling.data.tfl.gov.uk/usage-stats/320JourneyDataExtract01Jun2022-07Jun2022.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/321JourneyDataExtract08Jun2022-14Jun2022.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/322JourneyDataExtract15Jun2022-21Jun2022.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/323JourneyDataExtract22Jun2022-28Jun2022.csv',
 'https://cycling.data.tfl.gov.uk/usage-stats/324JourneyDataExtract29Jun2022-05Jul2022.csv']


for journey_data_url in tqdm(selected_journey_data):
    fname = journey_data_url.split('/')[-1]
    path = DATA_DIR + '/' + fname
    
    if os.path.exists(path):
        continue
        
    r = requests.get(journey_data_url)
        
    with open(path, 'a') as f:
        f.write(r.text)

100%|██████████████████████████████████████████| 5/5 [00:00<00:00, 24995.85it/s]


FOI_TFL_CYCLE_STATION_DATA_URL = 'https://www.whatdotheyknow.com/request/766040/response/1846134/attach/4/FOI%200689%202122.csv.txt'

foir_station_df = pd.read_csv(FOI_TFL_CYCLE_STATION_DATA_URL)


# {`Original Column Name`: `Desired Column Name in Output`}
column_name_mappings = {
    'rentalid': 'rental_id',
    'startdate': 'start_date',
    'enddate': 'end_date',
    'bikeid': 'bike_id',
    'duration': 'duration_s',
    'duration_seconds': 'duration_s',
    'startstationid': 'start_station_id',
    'endstationid': 'end_station_id',
    'endstationname': 'end_station_name',
    'startstationname': 'start_station_name',
}


# To validate station names and identifiers, they are logged in sets.
# They get pickled in the procedure.

if not os.path.exists(INTERIM_OUTPUT_DIR + 'unique_station_identifiers.pickle'):
    unique_station_identifiers = set()
else:
    with open(INTERIM_OUTPUT_DIR + 'unique_station_identifiers.pickle', 'rb') as f:
        unique_station_identifiers = pickle.load(f)

for csv_path in tqdm(os.listdir(DATA_DIR)):
    feather_path = csv_path.replace('.csv', '.feather')

    # Check whether we've already processed this csv file:
    if feather_path in os.listdir(DATA_DIR):
        continue

    if '.csv' in csv_path:
        journey_extract_df = pd.read_csv(DATA_DIR + '/' + csv_path)

        # Sometimes TfL might write StartStation ID, and in other files they write Start Station Id.
        # Let's just remove all spaces and convert to lower to try and by as agnostic as possible to future changes.
        journey_extract_df.columns = [col.replace(' ', '').strip().lower() for col in journey_extract_df.columns]

        journey_extract_df.rename(column_name_mappings, axis='columns', inplace=True)

        # There are a few essential columns that are required for every journey.
        essential_numeric_columns = ['bike_id', 'end_station_id', 'start_station_id', 'duration_s']

        # If all the essential columns are not available, ignore the csv file.
        if not all(col in journey_extract_df.columns for col in essential_numeric_columns):
            continue

        journey_extract_df.dropna(subset=essential_numeric_columns, inplace=True)
        journey_extract_df.reset_index(drop=True, inplace=True)

        # Reduce memory use by storing smaller integers as smaller types:
        for colname in essential_numeric_columns:
            journey_extract_df[colname] = journey_extract_df[colname].astype('int16')

        # Convert dates to proper formats:
        journey_extract_df['start_date'] = pd.to_datetime(journey_extract_df['start_date'], dayfirst=True)
        journey_extract_df['end_date'] = pd.to_datetime(journey_extract_df['end_date'], dayfirst=True)

        # Remove any weird journeys with impossible dates
        journey_extract_df = journey_extract_df[journey_extract_df['start_date'] < journey_extract_df['end_date']]
        journey_extract_df['duration_s'] = journey_extract_df['duration_s'].abs()

        # Next, journeys are stored to the data folder.
        # Unnecessary / repeating columns are excluded (like station names).
        journey_extract_df[
            essential_numeric_columns +
            ['rental_id', 'end_date', 'start_date']
            ].to_feather(DATA_DIR + '/' + feather_path)

        # Validate the stations in the journey data against those in the
        # docking station dataset
        for station_id_name_pair in list(journey_extract_df.groupby(['start_station_id', 'start_station_name']).groups):
            unique_station_identifiers.add(station_id_name_pair)
        for station_id_name_pair in list(journey_extract_df.groupby(['end_station_id', 'end_station_name']).groups):
            unique_station_identifiers.add(station_id_name_pair)

        # Pickle the sets to serve as a backup if the script is interrupted.
        with open(INTERIM_OUTPUT_DIR + 'unique_station_identifiers.pickle', 'wb') as f:
            pickle.dump(unique_station_identifiers, f)

100%|████████████████████████████████████████| 12/12 [00:00<00:00, 18957.31it/s]


journeys_station_df = pd.DataFrame(
    unique_station_identifiers,
    columns = ['station_id', 'station_name']
)

journeys_station_df = journeys_station_df.join(
    journeys_station_df.station_id.value_counts().rename('occurrence_count'),
    on='station_id'
)


journeys_station_df['station_first_name'] = journeys_station_df.station_name.str.split(',').str[0].str.replace('\d+', '', regex=True).str.strip().str.lower()
journeys_station_df['station_area'] = journeys_station_df.station_name.str.split(',').str[1].str.replace('\d+', '', regex=True).str.strip().str.lower()


journeydata_problematic_station_ids = []

for station_id, group in journeys_station_df[
    journeys_station_df.occurrence_count > 1
  ].groupby('station_id'):

  if not (
      len(group.station_first_name.unique()) == 1 
      or len(group.station_area.unique()) == 1
  ):
    journeydata_problematic_station_ids.append(station_id)


len(journeydata_problematic_station_ids)

0


journeys_station_df[
    journeys_station_df.station_id.isin(journeydata_problematic_station_ids)
].sort_values('station_id')


# Rename the columns to align with Python convention. 
# Add an `foir_` prefix to feature names that exist in the `Cycle Journey Extract` dataset.
foir_station_df.columns = ['go_live', 'foir_docking_station_name', 'foir_docking_station_id', 
                           'n_docking_points', 'latitude', 'longitude']

foir_station_df['latitude'] = foir_station_df['latitude'].astype(float)
foir_station_df['longitude'] = foir_station_df['longitude'].astype(float)
foir_station_df['foir_docking_station_id'] = foir_station_df['foir_docking_station_id'].astype('int16')
foir_station_df['n_docking_points'] = foir_station_df['n_docking_points'].astype('int16')


foir_station_df


# Count the occurrence of every station_id and include it as a column:
foir_station_df = foir_station_df.join(
    foir_station_df.foir_docking_station_id.value_counts().rename('id_occurrence_count'),
    on='foir_docking_station_id'
)

foir_problematic_station_ids = []

# Iterate through groups of docking stations with the same `foir_docking_station_id`
for foir_docking_station_id, group in foir_station_df[
    foir_station_df.id_occurrence_count > 1
  ].groupby('foir_docking_station_id'):
    
    # If they do not all have roughly the same coordinates, mark the station id as problematic. 
    if not(len(group.latitude.round(3).unique()) == 1 and len(group.longitude.round(3).unique()) == 1):
        foir_problematic_station_ids.append(foir_docking_station_id)


foir_station_df[
    foir_station_df.foir_docking_station_id.isin(foir_problematic_station_ids)
].sort_values('foir_docking_station_id')


full_station_df = foir_station_df.merge(
    journeys_station_df, 
    left_on='foir_docking_station_id', 
    right_on='station_id', 
    how='outer'
)


full_station_df['foir_station_first_name'] = (
    full_station_df.foir_docking_station_name.str.split(',').str[0].str.replace('\d+', '', regex=True).str.strip().str.lower()
)
full_station_df['foir_station_area'] = (
    full_station_df.foir_docking_station_name.str.split(',').str[1].str.replace('\d+', '', regex=True).str.strip().str.lower()
)


non_matching_ids = full_station_df[(
  (full_station_df.foir_station_first_name != full_station_df.station_first_name) &
  (full_station_df.foir_station_area != full_station_df.station_area)
)].foir_docking_station_id.unique()


all_problematic_station_ids = (
    set(journeydata_problematic_station_ids) 
    or set(foir_problematic_station_ids) 
    or set(non_matching_ids)
)


print(f"Exclude {len(all_problematic_station_ids)} docking station ids from the analysis.")

Exclude 6 docking station ids from the analysis.


foir_station_df = foir_station_df[
    (~foir_station_df.foir_docking_station_id.isin(all_problematic_station_ids)) 
]

foir_station_df = foir_station_df.sort_values(
    'go_live', 
    ascending=False
).groupby('foir_docking_station_id').first()

foir_station_df.drop(
    ['go_live', 'id_occurrence_count'], 
    axis='columns', 
    inplace=True
)

# Remove the foir column name prefix that was added earlier.
foir_station_df.columns = [col.replace('foir_', '') for col in foir_station_df.columns]


def plot_docking_stations(docking_station_df, zoom=8.5, title=None):
    return px.scatter_mapbox(
        docking_station_df.reset_index(),
        lat='latitude',
        lon='longitude',
        size='n_docking_points',
        color='n_docking_points',
        opacity=0.7,
        zoom=zoom,
        hover_data=['docking_station_name', 'foir_docking_station_id'],
        title=title
    )
plot_docking_stations(foir_station_df)


foir_station_df.loc[502, ['latitude', 'longitude']] = foir_station_df.loc[494, ['latitude', 'longitude']]


plot_docking_stations(foir_station_df)


foir_station_df.reset_index().to_feather(INTERIM_OUTPUT_DIR + 'FOI Station Data Cleaned.feather')


foir_station_df = pd.read_feather(INTERIM_OUTPUT_DIR + 'FOI Station Data Cleaned.feather').set_index('foir_docking_station_id')


journey_df = pd.concat([
    pd.read_feather(DATA_DIR + '/' + fname)
    for fname in os.listdir(DATA_DIR)
    if '.feather' in fname
])


print(f"{len(journey_df)} cycle journeys and {len(foir_station_df)} docking stations have been loaded.")

1475810 cycle journeys and 797 docking stations have been loaded.


# The station data frame contains a row for every cycle hire station point on the TfL cycle hire network.
# Join this onto the journey data frame, and count the most popular journeys taken on the network.

def get_journey_counts_from_journey_data(df):
    """
    Count the unique journeys in a dataframe containing journey row data, and return the output.
    :param df: A DataFrame containing journey data.  Requires the columns `start_station_id` and `end_station_id`
    :return:
    """
    global foir_station_df
    
    journey_count_df = df[['start_station_id', 'end_station_id']].value_counts().to_frame()
    journey_count_df.columns = ['Journey Count']
    journey_count_df.reset_index(inplace=True)
    
    # Exclude journeys that start and end at the same station:
    journey_count_df = journey_count_df[journey_count_df['start_station_id'] != journey_count_df['end_station_id']]

    # Join onto the station data:
    journey_count_df = journey_count_df.join(
         foir_station_df,
         on='start_station_id'
    ).join(
         foir_station_df,
         on='end_station_id',
         rsuffix='_destination'
    ).dropna(
        subset=['docking_station_name', 'docking_station_name_destination'] 
    )
    
    return journey_count_df


journey_counts = get_journey_counts_from_journey_data(journey_df)


journey_count_sample = journey_counts[
    journey_counts['Journey Count'] > journey_counts['Journey Count'].quantile(0.9)
]


OSRM_BASE_URL = "https://router.project-osrm.org/route/v1"

def osrm_route_directions(
    origin_lat: float, origin_lon: float, 
    destination_lat: float, destination_lon: float,
    profile: str='bike'
    ):
    """
    Return a linestring bicycle route between an origin and destination using the Open Source Routing Machine.
    See: https://project-osrm.org/docs/v5.24.0/api/#

    :param origin_lat: The latitude for the starting point's coordinates (WGS84).
    :param origin_lon: The longitude for the starting point's coordinates (WGS84).
    :param destination_lat: The latitude for the destination point's coordinates (WGS84).
    :param destination_lon: The longitude for the destination point's coordinates (WGS84).
    :param profile: The OSRM profile to route with.  This can be `car`, `bike`, or `foot` (for pedestrian directions).
    :return: A LineString representing the recommended OSRM route between the origin and destination.
    """
    
    # Format the OSRM route URL and get the result from the open OSRM endpoint
    coord_string = f"{origin_lon},{origin_lat};{destination_lon},{destination_lat}"
    url = f"{OSRM_BASE_URL}/{profile}/{coord_string}"
    routes = requests.get(url).json()['routes']
    
    # OSRM can return multiple routes.  Take the first line, convert it to Shapely geometry and return it.
    if len(routes) > 0:
        geometry = polyline.decode(routes[0]['geometry'], geojson=True)
        return shapely_geometry.LineString(geometry)
    return None


if not os.path.exists(INTERIM_OUTPUT_DIR + 'station_pair_route_geometries.pickle'):
    station_pair_route_geometries = {}
else:
    with open(INTERIM_OUTPUT_DIR + 'station_pair_route_geometries.pickle', 'rb') as f:
        station_pair_route_geometries = pickle.load(f)


for ix, row in tqdm(journey_count_sample.reset_index(drop=True).iterrows()):
    if (row.start_station_id, row.end_station_id) in station_pair_route_geometries:
        continue
        
    station_pair_route_geometries[(row.start_station_id, row.end_station_id)] = osrm_route_directions(
        row['latitude'], row['longitude'], 
        row['latitude_destination'], row['longitude_destination']
    )
    # Pickle the sets to serve as a backup if the script is interrupted.
    if ix%100 == 0:
        with open(DATA_DIR + 'station_pair_route_geometries.pickle', 'wb') as f:
            pickle.dump(station_pair_route_geometries, f)

21929it [00:00, 51375.01it/s]


station_id_cross_with_geometry = gpd.GeoDataFrame.from_dict(
    station_pair_route_geometries, orient='index', columns=['geometry'], geometry='geometry', crs='EPSG:4326'
)

/opt/homebrew/Caskroom/miniforge/base/lib/python3.10/site-packages/pandas/core/dtypes/cast.py:128: ShapelyDeprecationWarning:

The array interface is deprecated and will no longer work in Shapely 2.0. Convert the '.coords' to a numpy array instead.


station_id_cross_with_geometry.plot()

<AxesSubplot: >


journey_df['is_weekend'] = (journey_df.start_date.dt.weekday > 4)


journey_df['duration_min'] = journey_df.duration_s / 60


def distribution_plot_weekend_vs_weekday(df, y, xlim):
    axes = plt.gca()

    for weekend_flag in [True, False]:
        label = 'During the Weekend' if weekend_flag else 'During Weekdays'
        axes = df[df.is_weekend == weekend_flag].plot(
            y=y, 
            kind='kde',
            xlim=xlim, 
            ax=axes,
            label=label
        )

    return axes


def plot_cycle_hire_duration_distribution():
    distribution_plot_weekend_vs_weekday(journey_df, 'duration_min', (0, 80))
    plt.title("Distribution of Cycle Hire Duration")
    plt.xlabel("Rental Duration in Minutes")
    plt.show()
    
plot_cycle_hire_duration_distribution()


print(f"The median journey is {journey_df.duration_min.median()} minutes long.")

The median journey is 15.0 minutes long.


prop_journeys_under_30_min = ((journey_df.duration_min < 30).sum() / len(journey_df)).round(2)
print(f"{100 * prop_journeys_under_30_min}% of journeys are less than 30 minutes.")

86.0% of journeys are less than 30 minutes.


# Route geometries are indexed by `(start_station_id, end_station_id)` tuples.
# Derive this tuple and store it as a column.  Route geometry will be joined on this column.
journey_df['station_pair'] = journey_df[['start_station_id', 'end_station_id']].apply(tuple, 1)

journey_df = journey_df.join(
    station_id_cross_with_geometry,
    on='station_pair',
    how='inner'
)


# Convert to a GeoDataFrame to allow for distance calculations:
journey_df = gpd.GeoDataFrame(journey_df, geometry='geometry', crs='EPSG:4326')
journey_df['distance_m'] = journey_df.to_crs('EPSG:27700').geometry.length

MILES_TO_METERS = 1609.34 # A coefficient to convert between distances in miles and meters. 
journey_df['distance_mi'] = journey_df['distance_m'] / MILES_TO_METERS


def plot_cycle_hire_distance_distribution():
    distribution_plot_weekend_vs_weekday(journey_df, 'distance_mi', (0, 6))
    plt.title("Distribution of Cycle Journey Distance*")
    plt.xlabel("Distance of the recommended route (mi)")
    plt.figtext(0.1, -0.1, "*Assuming bicycles take the OSRM-recommended route between two stations.")
    plt.show()
    
plot_cycle_hire_distance_distribution()


print(f"The median journey is {round(journey_df.distance_mi.median(), 2)} miles long")

The median journey is 1.56 miles long


MPS_TO_MPH = 2.23694 # A coefficient to convert between meters-per-second and miles-per-hour


journey_df['speed_mph'] = MPS_TO_MPH * journey_df['distance_m'] / journey_df['duration_s']


journey_df = journey_df[journey_df['speed_mph'] < 30]


def plot_cycle_journey_speed_distribution(): 
    distribution_plot_weekend_vs_weekday(journey_df, 'speed_mph', (0, 30))
    plt.title("Distribution of Average Cycle Journey Speed*")
    plt.xlabel("Average Cycle Journey Speed (mph)")
    plt.figtext(0.1, -0.1, "*Assuming bicycles take the OSRM-recommended route between two stations.")
    plt.show()
    
plot_cycle_journey_speed_distribution()


journey_df.is_weekend.value_counts()

False    729144
True     252459
Name: is_weekend, dtype: int64


journey_counts_df = get_journey_counts_from_journey_data(
    journey_df
)


journey_counts_df.index = journey_counts_df[['start_station_id', 'end_station_id']].apply(tuple, 1)
journey_counts_df = journey_counts_df.join(station_id_cross_with_geometry)


journey_counts_gdf = gpd.GeoDataFrame(journey_counts_df.reset_index(drop=True), geometry='geometry')
journey_counts_gdf.to_file(INTERIM_OUTPUT_DIR + 'TfL Cycle Journey Counts.geojson')


with open('./config/route_popularity_kepler_config.json') as f:
    route_popularity_kepler_config = json.load(f)
    
route_popularity_kepler_chart = keplergl.KeplerGl(
    height=500,
    data={
        'p5nn3waud': journey_counts_gdf[
            journey_counts_gdf['Journey Count'] > 50
        ]
    },
    config=route_popularity_kepler_config
)

User Guide: https://docs.kepler.gl/docs/keplergl-jupyter


route_popularity_kepler_chart.save_to_html(
    file_name=REPORT_OUTPUT_DIR + 'Route Popularity Visualisation.html',
    data={
        'p5nn3waud': journey_counts_gdf[
            journey_counts_gdf['Journey Count'] > 50
        ]
    },
    config=route_popularity_kepler_config,
    read_only=True
)

Map saved to ./Report Output/Route Popularity Visualisation.html!


def extract_journey_counts_by_station_by_period(df, start_or_end, time_group_freq='15min'):
    """
    Convert a dataframe of journeys, to one with journey counts, by hour and station id.
    Specify whether to count journeys `start`ing or `end`ing with the `start_or_end` parameter.
    :param df: A dataframe of journeys.  Requires the following columns: `start_or_end`_station_id, rental_id, `start_or_end`_time.
    :param start_or_end: Either `start` or `end`, depending on whether to count journeys started, or journeys ended.
    :param time_group_freq: A period over which to count the journeys. Journeys are counted by station_id, time_group cells.
    :return: A DataFrame containing hire counts for every station and period combination.
    """
    df['date'] = df.start_date.dt.date

    # Extract the time of the rental / return to the nearest 15 minutes
    df.loc[:, f'{start_or_end}_time'] = df[f'{start_or_end}_date'].dt.round(time_group_freq).dt.time
    
    # Count returns/rentals by day, time and station id
    journeys_by_station_by_day_by_period = df.groupby(
        ['date', f'{start_or_end}_time', f'{start_or_end}_station_id']
    )[['rental_id']].count()
    
    # Fill in days and times with no rentals with zero
    journeys_by_station_by_day_by_period_pivot = pd.pivot_table(
        journeys_by_station_by_day_by_period,
        values='rental_id',
        columns=[f'{start_or_end}_station_id'],
        index=['date', f'{start_or_end}_time']
    ).fillna(0)
    
    journeys_by_station_by_day_by_period_unpivot = journeys_by_station_by_day_by_period_pivot.reset_index().melt(
        id_vars=['date', f'{start_or_end}_time'],
        value_name=f'n_hires_{start_or_end}ed_per_period'
    )
    
    # Average by time period
    journeys_by_station_by_period = journeys_by_station_by_day_by_period_unpivot.reset_index().groupby(
        [f'{start_or_end}_time', f'{start_or_end}_station_id']
    ).mean(numeric_only=True)[[f'n_hires_{start_or_end}ed_per_period']]
    
    return journeys_by_station_by_period.reset_index()


weekday_journey_df = journey_df[~journey_df.is_weekend].reset_index(drop=True)
weekend_journey_df = journey_df[journey_df.is_weekend].reset_index(drop=True)


weekday_started_journeys_by_station_by_time = extract_journey_counts_by_station_by_period(weekday_journey_df, 'start')
weekend_started_journeys_by_station_by_time = extract_journey_counts_by_station_by_period(weekend_journey_df, 'start')


def plot_journeys_started_weekend_vs_weekday():
    weekday_started_journeys_by_station_by_time.groupby('start_time').mean()['n_hires_started_per_period'].plot(
        figsize=(12, 5),
        label='Weekday #hires started per period'
    )
    weekend_started_journeys_by_station_by_time.groupby('start_time').mean()[f'n_hires_started_per_period'].plot(
        figsize=(12, 5),
        label='Weekend #hires started per period',
    )
    plt.xticks([time(h*2, 0) for h in range(12)])
    plt.legend()
    plt.xlabel("Time")
    plt.title("Average number of TfL Cycle Hire Journeys Started per 15-minute Period per Docking Station")
    plt.ylabel("# Hires per Docking Station per 15-minute Period")
    plt.show()
    
plot_journeys_started_weekend_vs_weekday()


def get_journey_start_end_station_counts(df):
    """
    Count journeys for every start-/end-station pair.
    :param df: A DataFrame containing journeys.  Required columns: `start_station_id`, `end_station_id`, `start_date`, `duration_s`
    :return: A DataFrame with journey counts for every station pair.
    """

    # Ignore journeys that start and end at the same station.
    filtered_journey_df = df[
        (df['start_station_id'] != df['end_station_id'])
    ].dropna()

    filtered_journey_df['start_hour'] = filtered_journey_df['start_date'].dt.hour
    filtered_journey_df['end_hour'] = filtered_journey_df['end_date'].dt.hour
    filtered_journey_df['date'] = pd.to_datetime(filtered_journey_df['start_date'].dt.date)
    
    # Count journeys by their date and the hour within which they started.
    journeys_by_station_and_hour = filtered_journey_df.groupby(
        ['start_station_id', 'end_station_id', 'date', 'start_hour']
    )[['rental_id']].count()
    journeys_by_station_and_hour.reset_index(inplace=True)
    
    # Count the number of unique dates for which trips exist for every unique station pair:
    station_n_unique_dates = journeys_by_station_and_hour.groupby(
        ['start_station_id', 'end_station_id']
    )['date'].nunique()
    
    # Count up the number of trips within every station-pair-hour cell, over all dates:
    station_total_trips_by_hour = journeys_by_station_and_hour.groupby(
        ['start_station_id', 'end_station_id', 'start_hour']
    )[['rental_id']].sum()
    
    # Add in the average trip duration for every station pair.
    station_total_trips_by_hour = station_total_trips_by_hour.join(
        filtered_journey_df.groupby(['start_station_id', 'end_station_id', 'start_hour'])[['duration_s']].mean()
    )
    
    # Join the unique date counts, to allow calculations of an average number of hires per hour:
    station_total_trips_by_hour = station_total_trips_by_hour.join(
        station_n_unique_dates,
        how='inner'
    )
    station_total_trips_by_hour['Hires Started per hour'] = (
        station_total_trips_by_hour['rental_id'] / station_n_unique_dates.max()
    )
    
    # Add a convenience feature for duration in minutes:
    station_total_trips_by_hour['Duration (min)'] = station_total_trips_by_hour.duration_s / 60
    
    # Join on the station location data.
    journey_start_end_station_counts = station_total_trips_by_hour.reset_index().join(
        foir_station_df,
        on='start_station_id'
    ).join(
        foir_station_df[['longitude', 'latitude']],
        on='end_station_id',
        rsuffix='_end'
    )
    
    return journey_start_end_station_counts


journey_start_end_station_counts_weekend = get_journey_start_end_station_counts(journey_df[journey_df.is_weekend])
journey_start_end_station_counts_weekday = get_journey_start_end_station_counts(journey_df[~journey_df.is_weekend])


with open('./config/kepler_arc_config.json') as f:
    kepler_arc_plot_config = json.load(f)


weekend_journey_arc_plot = keplergl.KeplerGl(
    height=500,
    data={'pj760hy89': journey_start_end_station_counts_weekend},
    config=kepler_arc_plot_config
)

User Guide: https://docs.kepler.gl/docs/keplergl-jupyter


weekend_journey_arc_plot.save_to_html(
    file_name=REPORT_OUTPUT_DIR + 'Weekend Journeys Arc Visualisation.html',
    data={'pj760hy89': journey_start_end_station_counts_weekend},
    config=kepler_arc_plot_config,
    read_only=True
)

Map saved to ./Report Output/Weekend Journeys Arc Visualisation.html!


weekday_journey_arc_plot = keplergl.KeplerGl(
    height=500,
    data={'pj760hy89': journey_start_end_station_counts_weekday},
    config=kepler_arc_plot_config
)

User Guide: https://docs.kepler.gl/docs/keplergl-jupyter


weekend_journey_arc_plot.save_to_html(
    file_name=REPORT_OUTPUT_DIR + 'Weekday Journeys Arc Visualisation.html',
    data={'pj760hy89': journey_start_end_station_counts_weekday},
    config=kepler_arc_plot_config,
    read_only=True
)

Map saved to ./Report Output/Weekday Journeys Arc Visualisation.html!


weekday_journey_arc_plot

KeplerGl(config={'version': 'v1', 'config': {'visState': {'filters': [{'dataId': ['pj760hy89'], 'id': 'xtkquj8…


weekend_journey_arc_plot

KeplerGl(config={'version': 'v1', 'config': {'visState': {'filters': [{'dataId': ['pj760hy89'], 'id': 'xtkquj8…


journey_df['date'] = pd.to_datetime(journey_df['start_date'].dt.date)


journey_df.date.value_counts().head(1)

2022-06-21    36539
Name: date, dtype: int64


single_day_journey_df = journey_df[journey_df.date == '2022-06-21']


linestring_featureset = []

for ix, row in tqdm(single_day_journey_df.iterrows()):
    
    # Draw out the coordinates from the linestring.
    # Calculate the distances between every successive coordinate.
    def split_line_segments(line):
        return list(map(LineString, zip(line.coords[:-1], line.coords[1:])))
    
    distances = gpd.GeoSeries(
        split_line_segments(row.geometry),
        crs='EPSG:4326'
    ).to_crs('EPSG:27700').length
    
    # Set the first coordinate as zero distance travelled:
    distances = np.array([0] + distances.values.tolist())
    
    if distances.sum() == 0:
        continue
    
    # Use the calculated distances to interpolate a number of seconds passed.
    # Assume the bicycle travels at a constant speed throughout the trip.
    time_passed_s = (distances / distances.sum() * row.duration_s).cumsum()
    
    # Add the seconds passed to the start date of the trip.
    timestamps = (row.start_date + pd.to_timedelta(time_passed_s, unit='s')).map(pd.Timestamp.timestamp)
    
    # Compile the coordinates, together with the timestamps into LineString features.
    route = row.geometry.coords
    if not any(pd.isna(time_passed_s)):
        linestring_featureset.append(
            Feature(
                geometry = LineString(zip(
                    np.array(route)[:,0],
                    np.array(route)[:,1],
                    np.zeros(len(route)),
                    timestamps.astype(int).values.tolist()
                )),
                properties = {
                    'id': row.rental_id,
                    'bike_id': row.bike_id,
                }
            )
        )

36539it [07:14, 84.13it/s]


collection = FeatureCollection(linestring_featureset)
with open(INTERIM_OUTPUT_DIR + "single_day_trip_features.json", "w") as f:
    f.write('%s' % collection)


with open('./config/trip_visualisation_kepler_config.json') as f:
    trip_visualisation_kepler_config = json.load(f)


trip_plot = keplergl.KeplerGl(
    height=500,
    data={'5ankfi16j': dict(collection)},
    config=trip_visualisation_kepler_config
)

User Guide: https://docs.kepler.gl/docs/keplergl-jupyter


trip_plot.save_to_html(
    file_name=REPORT_OUTPUT_DIR + 'Single Day Trip Visualisation.html',
    data={'5ankfi16j': dict(collection)}, 
    config=trip_visualisation_kepler_config,
    read_only=True
)

Map saved to ./Report Output/Single Day Trip Visualisation.html!


REPORT_OUTPUT_DIR + 'Single Day Trip Visualisation.html'

'./Report Output/Single Day Trip Visualisation.html'


# Start stations translate to bicycle rentals.  End stations are returns.
docking_station_rentals = journey_df.groupby(
    ['date', 'start_station_id']
)[['rental_id']].count().reset_index()

docking_station_returns = journey_df.groupby(
    ['date', 'end_station_id']
)[['rental_id']].count().reset_index()

docking_station_rentals.columns = ['date', 'foir_docking_station_id', 'n_rentals']
docking_station_returns.columns = ['date', 'foir_docking_station_id', 'n_returns']

docking_station_returns.set_index(['date', 'foir_docking_station_id'], inplace=True)
docking_station_rentals.set_index(['date', 'foir_docking_station_id'], inplace=True)

# Join the rentals and returns onto the same dataframe
docking_station_movements = docking_station_returns.join(docking_station_rentals, how='outer')


docking_station_movements['daily_surplus'] = (
    docking_station_movements.n_returns.fillna(0) - docking_station_movements.n_rentals.fillna(0)
)


docking_station_movements['daily_surplus'].sum()

0.0


def plot_bicycle_surplus_distribution():
    axes = plt.gca()

    docking_station_movements.plot(
        y='daily_surplus', 
        kind='kde',
        xlim=(-30, 30), 
        ax=axes,
    )

    plt.title("The Distribution of Daily Bicycle Surplus by Docking Station")
    plt.xlabel("Bicycle Surplus (# of Bicycles)")
    plt.show()
    
plot_bicycle_surplus_distribution()


average_daily_surpluse = docking_station_movements.daily_surplus.clip(lower=0).mean()


print(f"On average, for every docking station on the network, {round(average_daily_surpluse, 2)} bicycles need to be manually moved to another station at the end of the day to balance out the network.")

On average, for every docking station on the network, 3.09 bicycles need to be manually moved to another station at the end of the day to balance out the network.


average_docking_station_movements = docking_station_movements.reset_index().groupby('foir_docking_station_id').mean(numeric_only=True)
average_docking_station_movements = average_docking_station_movements.join(foir_station_df)
average_docking_station_movements.dropna(subset=['n_docking_points'], inplace=True)


fig_average_daily_station_surplus = px.scatter_mapbox(
    average_docking_station_movements,
    lat='latitude',
    lon='longitude',
    color='daily_surplus',
    hover_data=['docking_station_name'],
    size='n_docking_points',
    size_max=15,
    color_continuous_scale='RdYlGn',
    mapbox_style='dark',
    range_color=[-20, 20],
    opacity=0.9,
    zoom=11,
    height=700,
    title='Average Daily Bicycle Surplus by Docking Station'
)
fig_average_daily_station_surplus.show()


len(journey_df)

981603


journey_df.head()


len(foir_station_df)

797


foir_station_df.head()


plot_docking_stations(foir_station_df, zoom=10.5, title='Cycle Hire Docking Stations')


len(station_id_cross_with_geometry)

74201


station_id_cross_with_geometry.plot()

<AxesSubplot: >


plot_cycle_hire_duration_distribution()


plot_cycle_hire_distance_distribution()


plot_cycle_journey_speed_distribution()


print(f"The median journey's average speed: {round(journey_df.speed_mph.median(), 2)} mph")

The median journey's average speed: 8.16 mph


plot_journeys_started_weekend_vs_weekday()


plot_bicycle_surplus_distribution()


print(f"On average, for every docking station on the network, {round(average_daily_surpluse, 2)} bicycles need to be manually moved to another station at the end of the day to balance out the network.")

On average, for every docking station on the network, 3.09 bicycles need to be manually moved to another station at the end of the day to balance out the network.


fig_average_daily_station_surplus.show()

	go_live	foir_docking_station_name	foir_docking_station_id	n_docking_points	latitude	longitude
0	Jul-10	River Street, Clerkenwell	1	19	51.529200	-0.109971
1	Jul-10	Phillimore Gardens, Kensington	2	37	51.499600	-0.197574
2	Jul-10	Christopher Street, Liverpool Street	3	32	51.521300	-0.084606
3	Jul-10	St. Chad's Street, King's Cross	4	23	51.530100	-0.120974
4	Jul-10	Sedding Street, Sloane Square	5	27	51.493100	-0.156876
...	...	...	...	...	...	...
832	Dec-20	Canada Water Station, Rotherhithe	844	35	51.498439	-0.049150
833	Jan-21	Clapham Common Station, Clapham Common	355	21	51.461388	-0.139221
834	Jan-21	Gauden Road, Clapham	808	28	51.464995	-0.130909
835	Jan-21	Bermondsey Station, Bermondsey	845	30	51.497995	-0.063471
836	Apr-21	London Fields, Hackney Central	614	21	51.541190	-0.058826

	go_live	foir_docking_station_name	foir_docking_station_id	n_docking_points	latitude	longitude	id_occurrence_count
342	Oct-10	Oval Way, Vauxhall	355	21	51.486600	-0.117286	2
833	Jan-21	Clapham Common Station, Clapham Common	355	21	51.461388	-0.139221	2
459	Mar-12	Castalia Square, Cubitt Town	474	39	51.498100	-0.011457	2
831	Dec-20	Gascoyne Road, Victoria Park	474	20	51.541515	-0.038558	2
467	Mar-12	Thornfield House, Poplar	482	26	51.509800	-0.023770	2
824	Jul-20	Exhibition Road Museums 2, South Kensington	482	14	51.499900	-0.174554	2
574	Oct-13	Bradmead, Nine Elms	614	40	51.478200	-0.144691	2
836	Apr-21	London Fields, Hackney Central	614	21	51.541190	-0.058826	2
775	Feb-16	Bevington Road, North Kensington	805	27	51.520074	-0.206338	2
822	Jun-20	The Metropolitan, Portobello	805	18	51.520509	-0.200805	2
768	Dec-15	Stockwell Roundabout, Stockwell	808	33	51.473500	-0.122556	2
834	Jan-21	Gauden Road, Clapham	808	28	51.464995	-0.130909	2

	bike_id	end_station_id	start_station_id	duration_s	rental_id	end_date	start_date	is_weekend	duration_min	station_pair	geometry	distance_m	distance_mi	speed_mph	date
0	18241	500	542	420	121447141	2022-06-16 14:14:00	2022-06-16 14:07:00	False	7.0	(542, 500)	LINESTRING (-0.03376 51.51405, -0.03406 51.514...	1996.817814	1.240768	10.635147	2022-06-16
11327	5642	500	542	420	121597697	2022-06-20 08:01:00	2022-06-20 07:54:00	False	7.0	(542, 500)	LINESTRING (-0.03376 51.51405, -0.03406 51.514...	1996.817814	1.240768	10.635147	2022-06-20
11328	8835	500	542	420	121455637	2022-06-16 17:25:00	2022-06-16 17:18:00	False	7.0	(542, 500)	LINESTRING (-0.03376 51.51405, -0.03406 51.514...	1996.817814	1.240768	10.635147	2022-06-16
26965	18328	500	542	420	121497723	2022-06-17 15:38:00	2022-06-17 15:31:00	False	7.0	(542, 500)	LINESTRING (-0.03376 51.51405, -0.03406 51.514...	1996.817814	1.240768	10.635147	2022-06-17
64930	14907	500	542	480	122241591	2022-07-04 09:22:00	2022-07-04 09:14:00	False	8.0	(542, 500)	LINESTRING (-0.03376 51.51405, -0.03406 51.514...	1996.817814	1.240768	9.305753	2022-07-04

	docking_station_name	n_docking_points	latitude	longitude
foir_docking_station_id
1	River Street, Clerkenwell	19	51.5292	-0.109971
2	Phillimore Gardens, Kensington	37	51.4996	-0.197574
3	Christopher Street, Liverpool Street	32	51.5213	-0.084606
4	St. Chad's Street, King's Cross	23	51.5301	-0.120974
5	Sedding Street, Sloane Square	27	51.4931	-0.156876

P4DS (ODL1) Assignment 3¶

Data Analysis Project¶

Notebook template design: Brandon Bennett (2020.11.03, revised 2021.03.02)¶

Bicycle Sharing in London¶

Student: Rayno Willem Mostert¶

Email: od22rwm@leeds.ac.uk¶

Contents:¶

1. Project Plan¶

- The Data¶

- Project Aim & Objectives¶

- System Design¶

2. Program Code¶

- Objective 1: An Automated Data Load, Validation and Enrichment Script¶

- Objective 2: Illustrate the System Utility¶

- Objective 3: Analyse the Self-Sufficiency of the Network¶

3. Project Outcome¶

- Overview of Results¶

- Conclusion¶

README¶

FAQ:¶

This document is really long. Where's the good stuff?¶

What are all these folders?¶

Running the code¶

Where's the data?¶

Project Plan¶

The Data¶

Features¶

What is a Journey?¶

Concerns¶

Links to data sources:¶

Project Aim and Objectives¶

Specific Objective(s)¶

Objective 1: An Automated Data Load, Validation and Enrichment Script¶

Objective 2: Illustrate the System Utility¶

Objective 3: Analyse the balance of the Cycle Hire Network.¶

System Design¶

Architecture¶

Processing Modules and Algorithms¶

Program Code¶

Organising the Directory Structure¶

Objective 1: An Automated Data Load, Validation and Enrichment Script¶

Data Compilation¶

TfL Open Cycle Hire Journey Data¶

FOI Station Data¶

Data Cleaning¶

Cleaning the FOI Station Data:¶

Data Enrichment¶

Most Popular Journeys¶

Objective 2: Illustrate the System Utility¶

Rental Duration¶

Journey Distance¶

Journey Speed¶

Most Popular Journeys¶

Temporal Analysis¶

Journey Visualisation¶

Add some Animation¶

Objective 3: Analyse the Self-Sufficiency of the Network¶

Project Outcome¶

Overview of Results¶

Objective 1: An Automated Data Load, Validation and Enrichment Script¶

Objective 2: Analyse the System Utility¶

High-Level Metrics¶

Cycle Hire Duration¶

Cycle Journey Distance¶

Cycle Hire Speed¶

Most Popular Routes¶

🔍 Interactive Chart: Route Popularity¶

Temporal Analysis¶

🔍 Interactive Chart: Weekday Journey Popularity¶

🔍 Interactive Chart: Weekend Journey Popularity¶

⭐️🔍🚲 Interactive Chart: Plotting Moving Trips for a Single Day¶

Objective 3: Network Balance¶

Conclusion¶

Achievements¶

Limitations¶

Future Work¶

1. Project Plan ¶

- The Data ¶

- Project Aim & Objectives ¶

- System Design ¶

2. Program Code ¶

- Objective 1: An Automated Data Load, Validation and Enrichment Script ¶

- Objective 2: Illustrate the System Utility ¶

- Objective 3: Analyse the Self-Sufficiency of the Network ¶

3. Project Outcome ¶

- Overview of Results ¶

- Conclusion ¶

🔍 Interactive Chart: Route Popularity ¶

🔍 Interactive Chart: Weekday Journey Popularity ¶

🔍 Interactive Chart: Weekend Journey Popularity ¶

⭐️🔍🚲 Interactive Chart: Plotting Moving Trips for a Single Day ¶